﻿using System.Collections.Generic;
using System.Diagnostics.Contracts;
using System.Runtime.CompilerServices;
using System.Threading;
using MathService.DifferencialEquations.Numerical;

namespace System
{
    using Function = Func<double, double>;

    /// <summary>Класс методов-расширений для функции</summary>
    public static class FuncExtentions
    {


        public static Func<TArg2, Func<TArg1, TResult>> Carring<TArg1, TArg2, TResult>(this Func<TArg1, TArg2, TResult> f) { return y => x => f(x, y); }

        public static Function Add(this Function f, double a) { return a.Equals(0) ? f : x => f(x) + a; }

        public static Function Add(this Function f1, Function f2) { return x => f1(x) + f2(x); }

        public static Function ArgumentReverse(this Function f) { return x => f(-x); }

        public static Function ArgumentShift(this Function f, double x0)
        {
            return x0.Equals(0) ? f : x => f(x - x0);
        }

        public static Function ArgumentCompression(this Function f, double k, double b = 0)
        {
            return k.Equals(1) && b.Equals(0) ? f : x => f(k * x + b);
        }

        public static Function Divade(this Function f, double a)
        {
            return a.Equals(1) ? f : x => f(x) / a;
        }

        public static Function Divade(this Function f1, Function f2) { return x => f1(x) / f2(x); }

        /// <summary>Функция от функци q(f(x))</summary>
        /// <param name="f">Внутренняя функция</param>
        /// <param name="q">Внешняя функция</param>
        /// <returns>ФУнкция q от функции f</returns>
        public static Function Func(this Function f, Function q) { return x => q(f(x)); }

        /// <summary>Функция от функции f(q(x))</summary>
        /// <param name="f">Внешняя функция</param>
        /// <param name="q">Внутренняя функция</param>
        /// <returns>Функция f от функции q</returns>
        public static Function FuncFor(this Function f, Function q) { return x => f(q(x)); }

        public static Function GetAKF(this Function f, double deltaX, double x0 = 0.0)
        {
            return f.GetConvolution(f.ArgumentReverse(), deltaX, x0);
        }

        public static Function GetConvolution(this Function f, Function g, double deltaX, double x0 = 0.0)
        {
            var lv_Xmin = x0 - deltaX / 2;
            var lv_Xmax = x0 + deltaX / 2;
            return x => ((Function)(t => f(t) * g(x - t))).GetIntegralValue_Adaptive(lv_Xmin, lv_Xmax);
        }

        #region Численное дифференцирование

        public static Function GetDifferencial(this Function f, double dx = 0.0001, int n = 0) { return x => f.GetDifferencialValue(x, dx, n); }

        public static double GetDifferencialValue(this Function f, double x, double dx = 0.0001, int n = 0)
        {
            Contract.Requires<ArgumentOutOfRangeException>(n >= 0 && n <= 4,
                                                           "значение n должно быть в пределах 0 <= n <= 4");

            //if(n < 0)
            //    throw new ArgumentOutOfRangeException("n", "значение n должно быть в пределах 0 < n < 5");

            if(n == 0)
            {
                var lv_DeltaX = dx / 2.0;
                var lv_Xmin = x - lv_DeltaX;
                var lv_Xmax = x + lv_DeltaX;
                var lv_DeltaY = f(lv_Xmax) - f(lv_Xmin);
                return !dx.Equals(0.0)
                                   ? lv_DeltaY / dx
                                   : (lv_DeltaY > .0
                                                  ? double.PositiveInfinity
                                                  : (lv_DeltaY < .0 ? double.NegativeInfinity : double.NaN));
            }
            n--;
            var lv_Result = .0;
            for(var i = 0; i < 6; i++)
                lv_Result += Solover.Differential.diff_a[n, i] * f(x + (i * dx));
            return lv_Result / Solover.Differential.diff_b[n] / dx;
        }

        #endregion


        #region Численное интегрирование

        /// <summary>Определённый численный интеграл функции (метод трапеций)</summary>
        /// <param name="f">Интегрируемая функция</param>
        /// <param name="x1">Начальное значение интервала интегрирования</param>
        /// <param name="x2">КОнечное значение интервала интегрирования</param>
        /// <param name="dx">Шаг интегрирования</param>
        /// <returns>Значение результата численного интегрирования функции методом трапеций</returns>
        public static double GetIntegralValue(this Function f, double x1, double x2, double dx = 0.0001)
        {
            Contract.Requires(f != null);
            Contract.Requires(!dx.Equals(0));
            if(x1.Equals(x2)) return 0;

            var lv_Result = .0;
            var lv_DeltaX = x2 - x1;
            var N = Math.Abs(lv_DeltaX) / dx;
            var x = x1;
            var y = f(x);
            for(var i = 0; i < N; i++)
                lv_Result += y + (y = f(x += dx));
            return .5 * dx * lv_Result;
        }

        public static double GetIntegralValue(this Function f, Func<double, double, double> Core, double x1,
                                              double x2, double dx = 0.0001)
        {
            Contract.Requires(f != null);
            Contract.Requires(Core != null);
            Contract.Requires(!dx.Equals(0));
            if(x1.Equals(x2)) return 0;

            var lv_DeltaX = x2 - x1;
            var N = Math.Abs(lv_DeltaX) / dx;
            var x = x1;
            var lv_Core = Core(f(x), x);
            var lv_Result = .0;
            for(var i = 0; i < N; i++)
                lv_Result += lv_Core + (lv_Core = Core(f(x += dx), x));
            return .5 * dx * lv_Result;
        }

        /// <summary>Численный расчёт определённого интеграла методом симпсона</summary>
        /// <param name="f">Интегрируемая функция</param>
        /// <param name="x1">Нижний предел интегрирования</param>
        /// <param name="x2">Верхний предел интегрирования</param>
        /// <param name="N">Число интервало интегрирования N > 2</param>
        /// <returns>Интеграл функции на отрезке метдом Симпсона</returns>
        public static double GetIntegralValue_Simpson(this Function f, double x1, double x2, int N = 100)
        {
            Contract.Requires(f != null);
            Contract.Requires(N > 2);
            if(x1.Equals(x2)) return 0;

            var dx = Math.Abs(x2 - x1) / N;
            var dx05 = dx / 2;
            var x = x1;
            var sum1 = f(x + dx05);
            var sum2 = .0;
            for(var i = 1; i < N; i++)
            {
                x += dx;
                sum1 += f(x + dx05);
                sum2 += f(x);
            }
            return dx / 6 * (f(x1) + 4 * sum1 + 2 * sum2 + f(x2));
        }

        /// <summary>Численный расчёт определённого интеграла методом адаптивного разбиения</summary>
        /// <param name="f">Интегрируемая функция</param>
        /// <param name="x1">Нижний предел интегрирования</param>
        /// <param name="x2">Верхний предел интегрирования</param>
        /// <param name="N">Начальное разбиение отрезка</param>
        /// <param name="eps">Точность вычисления интеграла</param>
        /// <returns>Адаптивный интеграл функции</returns>
        public static double GetIntegralValue_Adaptive(this Function f, double x1, double x2, int N = 2,
                                                       double eps = 0.000001)
        {
            Contract.Requires(f != null);
            Contract.Requires(N > 2);
            Contract.Requires(eps >= 0);
            if(x1.Equals(x2)) return 0;
            if(eps.Equals(0)) eps = double.Epsilon;

            var I1 = f.GetIntegralValue_Simpson(x1, x2, N);
            N *= 2;
            var I2 = f.GetIntegralValue_Simpson(x1, x2, N);
            return Math.Abs(I1 - I2) < eps
                               ? I2
                               : f.GetIntegralValue_Adaptive(x1, .5 * (x1 + x2), N, eps)
                                 + f.GetIntegralValue_Adaptive(.5 * (x1 + x2), x2, N, eps);
        }

        //public static double GetIntegralValue_InfiniteIntervals(this Function f, double x1, double x2, int N = 2, double eps = 0.000001)
        //{
        //    if(x2 < x1)
        //        return -f.GetIntegralValue_InfiniteIntervals(x2, x1, N, eps);

        //    var x1_inf = double.IsNegativeInfinity(x1);
        //    var x2_inf = double.IsPositiveInfinity(x2);

        //    if(x1_inf && x2_inf)
        //    {
        //        Func<double, double> F = t =>
        //            {
        //                var t2 = t * t;
        //                var T = 1 - t2;
        //                return f(t / T) * (1 + t2) / T / T;
        //            };
        //        return F.GetIntegralValue_Adaptive(0, 1, N, eps);
        //    }
        //    if(x1_inf)
        //    {
        //        Func<double, double> F = t =>
        //            {
        //                var t2 = t * t;
        //                return f(x2 - (1 - t) / t) / t2;
        //            };
        //        return F.GetIntegralValue_Adaptive(0, 1, N, eps);
        //    }
        //    if(x2_inf)
        //    {
        //        Func<double, double> F = t =>
        //            {
        //                var T = 1 - t;
        //                return f(x1 + t / T) / T / T;
        //            };
        //        return F.GetIntegralValue_Adaptive(0, 1, N, eps);
        //    }
        //    return f.GetIntegralValue_Adaptive(x1, x2, N, eps);
        //}


        /// <summary>Значение интеграла функции методом сплайнов</summary>
        /// <param name="f">Интегрируемая функция</param>
        /// <param name="x1">Начало отрезка интегрирования</param>
        /// <param name="x2">Конец отрезка интегрирования</param>
        /// <param name="N">Разбиение интервала интегрирования</param>
        /// <returns>Значение интеграла функции</returns>
        [Copyright("MachineLearning.ru", url = "http://www.machinelearning.ru/wiki/index.php?title=Применение_сплайнов_для_численного_интегрирования")]
        public static double GetIntegralValue_Spline(this Function f, double x1, double x2, int N = 1000)
        {
            Contract.Requires(f != null);
            Contract.Requires(N > 2);


            if(x1.Equals(x2)) return 0;

            var N1 = N - 1;
            var dx = (x2 - x1) / N;
            var f_data = new double[N + 1];
            double x;
            int i;
            for(i = 0; i <= N; i++)
            {
                x = x1 + dx * i;
                f_data[i] = f(x);
            }

            var alpha = new double[N1];
            alpha[0] = -1 / 4;
            var beta = new double[N1];
            beta[0] = f_data[2] - 2 * f_data[1] + f_data[0];
            //метод прогонки, прямой ход
            for(i = 1; i < N1; i++)
            {
                alpha[i] = -1 / (alpha[i - 1] + 4);
                beta[i] = (f_data[i + 2] - 2 * f_data[i + 1] + f_data[i] - beta[i - 1]) / (alpha[i - 1] + 4);
            }

            var c = new double[N1];
            c[N1 - 1] = (f_data[N] - 2 * f_data[N - 1] + f_data[N - 2] - beta[N1 - 1]) / (4 + alpha[N1 - 1]);
            //обратный ход
            for(i = N1 - 2; i >= 0; i--)
                c[i] = alpha[i + 1] * c[i + 1] + beta[i + 1];

            for(i = 0; i < N1; i++) c[i] = c[i] * 3 / (dx * dx);

            //считаем приближенное значение интеграла по формуле (9):
            x = (5 * f_data[0] + 13 * f_data[1] + 13 * f_data[N - 1] + 5 * f_data[N]) / 12;
            var tmp = 0.0;
            for(i = 2; i < N - 1; i++) tmp = tmp + f_data[i];
            x = (x + tmp) * dx - (c[0] + c[N1 - 1]) * dx * dx * dx / 36;
            return x;

            //считаем приближенное значение интеграла по формуле (6):
            //tmp = 0;
            //x = (f[0] + f[N]) / 2;
            //for(i = 1; i < N; i++) tmp = tmp + f[i];
            //x = (x + tmp) * tau;
            //tmp = 0;
            //for(i = 0; i < n; i++) tmp = tmp + c[i];
            //tmp = tmp * tau * tau * tau / 6;
            //return x;

        }

        public static Function GetIntegral(this Function f, double x0 = 0, double C = 0, double eps = 1e-6)
        {
            var integrator = new Integrator(f, x0, C);
            return x => integrator.GetValue(x, eps);
        }

        #endregion


        public static Function GetPeriodic(this Function f, double T, double x0 = 0.0)
        {
            if(!x0.Equals(0.0))
                f = f.ArgumentShift(x0);
            return x => f((x % T) + ((x < 0) ? T : 0));
        }

        public static double GetPower(this Function f, double x1, double x2)
        {
            return ((Function)(x =>
                         {
                             var num = f(x);
                             return num * num;
                         })).GetIntegralValue_Adaptive(x1, x2) / (x2 - x1);
        }

        public static IEnumerable<T> Sampling<T>(this Func<double, T> f, double x1, double x2, double dx)
        {
            var x = x1;
            do { yield return f(x += dx); } while(x <= x2);
        }

        public static TResult[] GetValues<TArgument, TResult>(this Func<TArgument, TResult> f, TArgument[] Argument)
        {
            return Argument.Function(f);
        }

        public static TResult[] GetValues<TResult>(this Func<double, TResult> f, double from, double to, double step)
        {
            var result = new List<TResult>((int)((to - from) / step));

            for(var i = 0; i < result.Capacity; i++)
                result.Add(f(from + i * step));

            return result.ToArray();
        }

        public static Function Inverse(this Function f) { return x => 1 / f(x); }

        public static Function Multiply(this Function f1, Function f2) { return x => f1(x) * f2(x); }

        public static Function Multiply(this Function f, double a) { return a.Equals(1) ? f : x => (f(x) * a); }

        public static Function Power(this Function f, double a)
        {
            return a.Equals(1) ? f : x => Math.Pow(f(x), a);
        }

        /// <summary>-f(x)</summary>
        public static Function Reverse(this Function f) { return x => -f(x); }

        public static Function SetParameter(this Func<double, double, double> f, double a, bool IsFirst = false)
        {
            return IsFirst ? (Function)(x => f(a, x)) : x => f(x, a);
        }

        public static Function Substract(this Function f1, Function f2) { return x => f1(x) - f2(x); }

        public static Function Substract(this Function f, double a) { return a.Equals(0) ? f : x => f(x) - a; }

        #region Nested type: FuncCalculator

        private class FuncCalculator<TArgument, TResult>
        {
            public struct func
            {
                public TArgument x;
                public Func<TArgument, TResult> f;
            }

            public readonly List<TArgument> Argument;
            public readonly List<TResult> Result;
            private int f_Complited;
            private int f_Count;

            public FuncCalculator()
            {
                Contract.Ensures(Argument != null);
                Contract.Ensures(Result != null);
                Argument = new List<TArgument>();
                Result = new List<TResult>();
            }

            public event EventHandler<EventArgs<List<TArgument>, List<TResult>>> OnComplited;

            [MethodImpl(MethodImplOptions.Synchronized)]
            public void Calculate(object State)
            {
                Contract.Requires(State != null);
                var lv_func = (func)State;
                Argument.Add(lv_func.x);
                Result.Add(lv_func.f(lv_func.x));
                f_Complited++;
                Check();
            }

            private void Check()
            {
                if((f_Count == f_Complited) && (OnComplited != null))
                    OnComplited(this, new EventArgs<List<TArgument>, List<TResult>>(Argument, Result));
            }

            [MethodImpl(MethodImplOptions.Synchronized)]
            public void SetCount(int Count)
            {
                f_Count = Count;
                Check();
            }
        }

        ///<summary>Вычислить значения функции параллельно</summary>
        ///<param name="f">Вычисляемая функция</param>
        ///<param name="Arguments">Область определения</param>
        ///<typeparam name="TArgument">Тип аргумента</typeparam>
        ///<typeparam name="TResult">Тип значения</typeparam>
        ///<returns>Массив значений функции</returns>
        public static TResult[] GetValuesParralel<TArgument, TResult>(this Func<TArgument, TResult> f,
                                                                      IEnumerable<TArgument> Arguments)
        {
            Contract.Requires(Arguments != null);
            var reset = new AutoResetEvent(false);
            var calculator = new FuncCalculator<TArgument, TResult>();
            calculator.OnComplited += (Sender, Args) => reset.Set();
            var count = 0;
            Arguments.Foreach(x =>
                                  {
                                      count++;
                                      ThreadPool.QueueUserWorkItem(calculator.Calculate,
                                                                   new FuncCalculator<TArgument, TResult>.func { x = x, f = f });
                                  });
            calculator.SetCount(count);
            reset.WaitOne();
            return calculator.Result.ToArray();
        }

        #endregion

        public class Integrator
        {
            private readonly object f_LockObject = new object();

            public double C { get; private set; }
            public double x0 { get; private set; }

            public Function f { get; private set; }

            public Integrator(Function f, double x0, double C = 0)
            {
                this.f = f;
                this.x0 = x0;
                this.C = C;
            }

            public double GetValue(double x, double eps = 1e-6)
            {
                lock(f_LockObject)
                {
                    if(x.Equals(x0))
                        return C;
                    var result = C += (x > x0
                                    ? f.GetIntegralValue_Adaptive(x0, x, eps: eps)
                                    : -f.GetIntegralValue_Adaptive(x, x0, eps: eps));
                    x0 = x;
                    return result;
                }
            }
        }
    }
}